#include "solver_omp_by_sav.h"

#include <math.h>
#include <stdlib.h>
#include <memory.h>
#include <conio.h>
#include <omp.h>
#include <Windows.h>
#include <stdio.h>
double* solver_omp_by_sav(double* u0, Area area, double epsilon, TwoDimFunc f)
{
	long workTime = GetTickCount();

	int n = area.n;
	double a = area.a, b = area.b;
	double* u_cur = new double[n * n];
	double* u_prev = new double[n * n];
	double h = (b - a) / n;
	double error;

	double* f_cur = new double[n * n];

	int numThreads=omp_get_max_threads();
	int threadID;

#pragma omp parallel shared (f_cur,a,h,numThreads) private(threadID)
	{
		threadID=omp_get_thread_num();
		for (int i = threadID; i < n; i+=numThreads)
		{
			for (int j = 0; j < n; ++j)
			{
				f_cur[i * n + j] = f(a + i * h, a + j * h);
			}
		}
	}

	memcpy(u_prev, u0, n * n * sizeof(double));
	memcpy(u_cur, u0, n * n * sizeof(double));

	do
	{ 
		numThreads=omp_get_max_threads();
#pragma omp parallel shared (u_cur, u_prev,h,n,f_cur,numThreads)  private(threadID)
		{
			threadID=omp_get_thread_num();
			for (int i = threadID + 1; i < n - 1; i+=numThreads)
			{
				for (int j = 1; j < n - 1; ++j)
				{
					u_cur[i * n + j] = 0.25 * 
						(u_prev[(i - 1) * n + j] + u_prev[(i + 1) * n + j] +
						u_prev[i * n + j - 1] + u_prev[i * n + j + 1] - 
						h * h * f_cur[i * n + j]);
				}
			}
		}

		double e;
		double* thredError = new double[numThreads];

#pragma omp parallel shared(thredError, u_cur, u_prev,n,numThreads) private(e)
		{	
			threadID=omp_get_thread_num();
			thredError[threadID] = abs(u_cur[threadID * n] - u_prev[threadID * n]);

			for (int i = threadID; i < n; i+=numThreads)
			{
				for (int j = 0; j < n; ++j)
				{
					e = abs(u_cur[i*n +j] - u_prev[i*n +j]);
					if (e > thredError[threadID]) thredError[threadID] = e;
				}
			}
		}

		error=0;
		for (int i=0;i<numThreads;i++)
		{
			if (error < thredError[i]) error = thredError[i];
		}

		double* tmp;
		tmp = u_cur; u_cur = u_prev; u_prev = tmp;

	} while (error > epsilon);

	delete[] u_cur;

	workTime= GetTickCount() - workTime;

	printf("Solver_omp_by_sav work time = %d\n",workTime);
	return u_prev;
}